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SEISMIC P-WAVE VELOCITY DERIVED FROM VIBRATOR CONTROL SYSTEM 

Field of the Invention 

This invention relates to seismic exploration, and more particularly to the mapping 
of underground features for oil and gas exploration. 

5 Background of the Invention 

Reliably estimating near-surface P-wave (longitudinal or compression wave) 
velocity is one of the key problems in land seismic exploration. Conventionally, uphole 
surveys have been used to acquire this near-surface velocity information. However, in 
the presence of rapidly varying near-surface geology, this method of estimation is 
10 inherently uncertain due to the sparse uphole sampling grid. 

Other conventional indirect methods to determine the near surface layers velocity 

t 

model have included refraction and shallow reflection surveys. These indirect 
measurements for near surface properties determination are conducted separately from the 
surveys directed towards imaging deep layers. Such refraction and shallow reflection 
15 surveys normally require further constraining factors to validate their results, and data 
from upholes has therefore been used for this extra validation. 

Correspondingly, data from refraction and shallow reflection surveys has been 
used to interpolate between sparsely located upholes for improving the near surface 
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velocity model. The velocity of the near surface layer is estimated from the direct 
arrivals from these two techniques. Sometimes other kinds of data are used for 
interpolation, such as ground penetrating radar (GPR) or resistivity. All these techniques 
are good aids for interpolation between upholes, but they add an extra cost component to 
5 the survey. 

A good knowledge of near-surface velocity macro model is vital for hydrocarbon 
reservoir exploration and characterization that utilize seismic data. This model is crucial 
for statics and depthing. However, the complexities of near surface layers make its 
determination, using direct measurements via uphole surveys, economically prohibitive. 
10 The problem that geophysicists face is how to interpolate between sparse near-surface 
velocity well controls knowing the fact ttiat near-surface velocity varies laterally with 
lithology that does not equally vary in all directions. Therefore, techniques are needed to 
estimate near-surface velocity using available data m order to minimize the associated risk 
resulting from an incomplete knowledge of the near-surface velocity model. 

15 Summary of the Invention 

It is therefore an object of the present invention to provide a method for 
estimating near-surface P-wave velocity that avoids the above-described difficulties of the 
prior art. 

The above and other objects are achieved by the present invention which, in one 
20 embodiment, is directed to a method of estimating P-wave velocity in a near-surface 
region of a land area, comprising the steps of gathering control data for the near-surface 
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region, gathering vibrator dynamic data generated in the near-surface region in response 
to vibrator action on the land area, and estimating the P-wave velocity in response to both 
the control data and the vibrator dynamic data. 

In a preferred embodiment, the control data is seismic data for the near-surface 
region generated in response to a shock in each of a plurality of upholes drilled in the 
land area, and the vibrator dynamic data includes both ground stiffness data and ground 
viscosity data. 

This approach of the present invention resolves the complexity of near-surface 
velocity model determination in areas characterized by sparse uphole controls and rapidly 
varying velocity. Moreover, the proposed technique does not add any cost for additional 
data acquisition because the needed data is readily available with seismic data acquired 
using vibrators. In addition, the proposed method requires less time than any of the 
conventional near-surface macro model estimation techniques. 

These and other objects, features and advantages of the present invention will be 
apparent from the following detailed description of the preferred embodiments taken in 
conjunction with the following drawings, wherein like reference numerals denote like 
elements. 

Brief Description of the Drawings 

Fig. 1 is a schematic cross section of an area to be surveyed having three upholes. 
Fig. 2 illustrates a conventional vibrator assembly mounted on a truck. 
Fig. 3 illustrates a basic model of vibrator theory. 
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Fig. 4 illustrates an idealized relationship of Vp/Vj to Poisson's ratio. 

Fig. 5 is a crossplot of calculated Vp versus 20 meter iso-depth uphole velocities. 

Fig. 6 illustrates a geostatistical integration model. 

Fig. 7 illustrates the steps of VALVE model building in accordance with the 
5 present invention. 

Fig. 8 illustrates the vibrator data geometry of a first test of the present invention. 
Fig. 9 illustrates near-surface velocity model building using geostatistical 
interpolation and integration in accordance with the present invention. 

Fig. 10 is a crossplot of derived velocity estunates measured independently at 
10 different times. 

Fig. 11 is a comparison between VALVE and measured uphole velocities. 
Fig. 12 illustrates three seismic sections along a first line (30 kilometers). 
Fig. 13 illustrates three seismic sections along a second line (30 kilometers). 
Fig. 14 shows the results of the 2-Layer model with real upholes posted. 
15 Fig. 15 shows the results of the 2-Layer model and geology incorporation with 

pseudo upholes posted. The contour interval is -5ms. 
Fig. 16 shows a seismic section along line CS-5. 

Fig. 17 is a basemap exhibiting the processed seismic data area and locations of 
the seismic sections. The background map shows the estimated P-wave velocity from 
20 vibrators. 

Fig. 18 shows three seismic sections obtained from the three statics models after 
the application of automatic statics. 
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Detailed Description of the Preferred Embodiments 

The present invention combines data from two sources to provide an improved 
estimate of the P-wave (compression wave) velocity in the near-surface layer. The first 
type of data is from upholes, also known as shot holes, tiiat are dug in the surface in an 
5 array spanning the area to be surveyed. A shock wave is generated on the surface and 
received by detectors attached to the bore-hole wall at pre-specified intervals. The 
obtained time versus depth data are used to estimate the P-wave velocity. Fig. 1 
illustrates in schematic cross-section an area to be surveyed using data from three 
illustrated upholes, and also indicates some of the variations in underground formations at 

10 the near-surface. Methods of using the data from such upholes for P-wave velocity 
estimation are well known and will be described herein only to the extent required to 
explain the present invention. 

The second type of data is from the vibrators that are conventionally used to 
measure vibrator and earth dynamics such as ground stiffness and ground viscosity. 

15 Vibrator estimated ground parameters and their respective geographic locations are 

routinely recorded in conjunction with the normal seismic acquisition of data for quality 
control purposes. Fig. 2 illustrates a conventional vibrator assembly mounted on a truck, 
with a schematic indication of variables measured in response to the vibrations impressed 
by these assemblies. 

20 A brief description of the underlying theory will now be given. The basic model 

for the vibrator data comes from the Lysmer equations, which describe the dynamic 
response of a rigid circular footmg to vertical motion, as in the model shown in Fig. 3. 
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This model shows a rigid circular footing of radius ro and mass m coupled to the elastic 
half-space and put into oscillation by an external periodic force Q, i.e. the action of the 
vibrator. 

Ground stiffness Kg and ground viscosity Dg are well known from Lysmer's 
5 Relations, as follows: 

where G is the shear modulus and v is Poisson's ratio. The units of Kg 
are (N/m). 

10 Dg = M£b- (poy^' 

where p is the mass density. The units of Dg are (N-sec/m). 

Poisson's ratio which is a parameter in both equations, can be 

expressed in terms of Vp (compression wave propagation velocity) and (shear 

15 wave propagation velocity) as follows: 

. = 0.5 - (V^ 
1 - (Vs/Vp)^ 

Also, the shear modulus G can be expressed in terms of mass density p 
and Vj by 
20 V3 = (G/pY^^ 

Combining these equations yields 
Kg = SToPiWs' -(V//Vp^)) and 
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Because Kg and Dg depend on three interdependent parameters Vp, and p, 
relating each of these ground parameters independently to ttie P-wave velocity Vp is not 
possible without knowledge of a further parameter. 

This further parameter comes from the critical damping for vertical oscillation in 
viscously damped systems defined such that the free displacement reaches equilibrium 
without oscillation. It is given by 

D, = 2(Kgm)^^ 

The critical damping is expressed in terms of Dg and a dimensionless mass ratio 
such that 

Dg = 0.425 

is the dimensionless mass ratio introduced by Lysmer given by 

Bz = Izii * m 
4 pro' 

Knowing both Kg and Dg, the equations yield 

-0- = ^^. 

1-p 2.89Kgro' 

Substitution into the definition of Kg finally yields 
V, = 0.85Kgro/Dg 

The previous derivation shows that both Kg and Dg can be used to eliminate the 
dependence on density and Poisson's ratio to compute V^. An idealized relationship of 
Vp/Vj to Poisson's ratio is shown in Fig. 4. In general, Poisson's ratio for cohesionless 



soils ranges from 0.25 to 0.35 and for cohesive soils from 0.35 to 0.45. The 
corresponding Vp/Vj ratio will vary from 1.73 to 2.08 for cohesionless soils and from 
2.08 to 3.32 for cohesive soils. The median of these ratios is 2.3. This value is used as 
a reasonable approximation of the Vp/Vg ratio for the near surface materials that will be 
sensed by the vibrator, as not all of these materials will be cohesive or cohesionless, but 
rather a combination of the two. 

Therefore, in accordance with the present invention, obtained from Kg and Dg 
will be multiplied by 2.3 and then correlated with Vp, seeking a linear relation. 

Table 1 below shows the correlation coefficients between collocated upholes and 
vibrator velocity attribute measurements at different iso-depths from the surface in a 
specified region selected for a test of the present invention. 

TABLE 1 



Depth (m) 


10 


20 


30 


40 


50 


60 


Vp 


0.71 


0.71 


0.69 


0.66 


0.56 


0.52 



The correlation coefficients pertaining to this attribute are greater than 0.65 down 
to a depth of 40 meters and decrease with depth, which confirms predictions. This due 
to the fact that estimates derived from vibrator measurements are only influenced by a 
relatively thin section of the earth surface. There is an abrupt change in correlation 
coefficients occurring at a depth of 50 meters, indicating that the depth of influence cutoff 
is between 40 and 50 meters. 
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Fig. 5 shows a crossplot of estimated Vp and 20 meter iso-depdi uphole velocities. 
The exhibited linear relation confirms the predictions. The plotted Vp is obtained by 
multiplying by 2.3, for the reasons given above. This fits the data well where the 
slope of the best fitted line is approximately 1. 
5 As described below, the method in accordance with the present invention was 

subjected to two field tests in order to evaluate its performance relative to other methods 
of velocity estimation. 

In the first test, uphole data of each iso-depth down to 50 meters was integrated 
with the vibrator velocity attribute using linear relations with the correlation coefficients 
10 shown in Table 1. Before integration, the vibrator velocity attribute values with kriged in 
order to form an exhaustive data set. Then the integration was done using Collocated 
Cokriging. The second test is described later in this specification. 

Thus, the present invention is directed to a novel technique, called herein 
"Vibrator Attribute Leading Velocity Estimation" (VALVE), for improving near-surface 
15 velocity determination using vibrator baseplate estimates of ground parameters. The first 
step is to develop a hypothesis for deriving P-wave velocity attribute from vibrator 
baseplate data. Geostatistics in data interpolation and integration is then applied. This is 
followed by generating an integrated 3D near-surface velocity model using data from a 
2,450 square kilometer area for the test. Finally, the results of applying this model on 
20 seismic data stacks are illustrated. 

Vibrator Baseplate Measurements As described above, while vibrating at any 
location on the ground, the vibrator exerts a force that is opposed by a counter force 
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from the ground. Simply, as it pushes against the earth, the vibrator feels the earth 
response to the applied force through the movements of the baseplate. Therefore, 
knowing the dynamics of the vibrator, that is the reaction mass and baseplate 
accelerations, estimates can be obtained for the underlying earth properties. Modern 
vibrator control systems estimate the actual ground force generated by vibrators relative 
to the theoretical input signal through measurements made at different parts of the 
vibrator, and thus, produce estimates of ground parameters at each vibration point (VP). 

As discussed in detail above, the theory of vibrations on the earth surface is well 
studied in the field of civil engineering, with the vertical oscillation of a footing resting 
on the surface of an elastic half-space being the civil engineering analogue of the vibrator 
and earth's dynamics. This theoretical framework is used to derive a velocity attribute of 
the near-surface from vibrator ground parameters measurements. This attribute has the 
same spatial sampling rate as the source grid in the seismic survey, which is considerably 
finer than the uphole grid. 

Geostatistical Estimation Tools Geostatistics is routinely used to predict reservoir 
parameters based on what is usually a limited number of available wells. Besides of 
describing spatial and temporal patterns, it is a good tool for multi-scaled data 
integration. In accordance with the present invention, the spatially well-sampled 3D 
seismic data is efficiently integrated with sparsely but vertically well sampled reservoir 
properties using geostatistical techniques (Fig. 6). In this case, seismic attributes are 
integrated with direct measurements, made in the wells, to improve the reliability of 
reservoir properties estunates away from the wells. Consequently, this integration 
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approach can be applied at any level within the earth layers (for example to improve a 
near surface layers velocity model obtained from direct measurements such as upholes) 
provided that other related data components are available. 

Kriging is the basic geostatistical interpolation tool. It is a local estimation 
technique that provides unbiased estimates witii minimum variance. It is also known as 
BLUE (Best Linear Unbiased Estimator) that provides optimal interpolation. It is unlike 
traditional interpolation techniques that depend on data values, because it incorporates a 
model of spatial correlation which makes it reliably honor the geologic features. 

Improvements in estimating missing spatial or temporal points are obtained when 
the primary measurements are integrated with related secondary attributes. In practice, 
the primary measurements are sparsely sampled compared to the densely sampled 
secondary information. There are different forms of kriging that allow performing this 
task. Among these are kriging with strata, simple kriging with varying local means, 
kriging with an external drift and cokriging. The latter besides kriging will be used to 
build the subsequent near-surface velocity models in the embodiments of the present 
invention. 

VALVE Near-Surface 3D Velocitv Model Building 

Near-surface 3D velocity model building using VALVE in accordance with the 
present invention consists of five steps, generally shown in Fig. 7: 

Step 1) Data base building. Three data components are required to generate a 
VALVE model. These are: 

a) Adequate upholes to ensure good statistical representation of the study area 
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near-surface velocity variations. 

b) Vibrator baseplate measurements of ground parameters: stiffness and 

viscosity. 

c) Surface elevation information. 

It should be rioted that the third data component is required if such a model is 
intended to be used for calculating seismic statics from surface to Seismic Reference 
Datum (SRD). 

Step 2) Uphole data preparation and quality assessment. Upholes represent the 
source of primary P-wave velocity data in accordance with the invention. In the test of 
an embodiment as described herein, 463 upholes were available in the study area with 
variable penetration depths. 444 of these were used as part of the primary dataset used to 
construct the VALVE model. The remaining 19 served as a control data set to compare 
against the predicted layer velocities. 

A depth versus time plot was constructed for each uphole to assist in detecting 
anomalous samples. Theoretically, travel time should increase with depth. However, 
this criterion is sometimes not satisfied due to various reasons. For example, cavities or 
fractures introduced by the drilling process can result in anomalous travel times. Another 
source of error could be experimental such as a wrongly picked arrival time or an 
incorrect depth. Errors caused by these factors normally manifest themselves on the 
depth-time plots, and thus can be corrected provided that enough control points are 
available. 

In general, uphole surveys cannot provide perfect results. For instance, lateral 
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changes in near-surface geology can cause inaccurate measurements depending on the 
magnitude of the source or the receiver offset from the well. The assumption of a 
straight raypath from source to receiver can be another source of errors, especially with 
heterogeneous near-surface layers. Nevertheless, when the source offset from the 
borehole is relatively small, as was true in the test, the magnitude of these errors is 
minimal. 

Finally, average P-wave velocities in iso-depth markers from the surface are 
calculated from the uphole data. These uphole average velocity markers are used for 
subsequent correlation and integration with the derived velocity attribute. 

Step 3) Vibrator ground parameters preparation and quality assessment 

Vibrator estimated ground parameters and their respective geographic locations are 
usually readily available from seismic surveys that use vibrators as sources. The 
measurements used herein were obtained from a 3D seismic survey conducted in the 
study area. The spatial distribution of such data was the same as the distribution of 
vibration points (VP's) in the survey design (Fig. 8). The VP's were acquired along 720 
meters apart east-west lines and 480 meters apart north-south lines. The VP's along each 
line were spaced at 60 meter intervals. Five vibrators were used at each VP spaced by a 
12 meter interval following an east-west pattern. 

Before deriving a velocity attribute from the vibrator data, the quality of the latter 
must be assessed. Since vibrator data are gathered in time-order, a time series 
presentation may give some indication about the behavior of these data. If the common 
cause of variations of vibrator measurements is due to the reaction to surface conditions. 
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then the majority of the vibrators operating at the same tune would provide approximately 
similar measurements. On the other hand, any specific variation in this process is 
assumed to be mainly due to problems associated with particular vibrators. Problems 
occurring with vibrators could be mechanical (i.e. related to the vibrator itself) or 
5 instrumental (related to the control system). Data measured by a vibrator that has a 
problem will deviate from the majority pattern, and data belonging to any vibrator that 
showed obvious deviation from the majority were removed from the analyzed data set. 

Acceptable measured data occur within three standard deviations from both sides 
of the mean. Therefore, pursuant to the previous visual quality assessment step, data 

10 samples that fall outside the three-sigma limits were removed from the data set. This 
operation was performed also on a daily basis. Following this step, the P-wave velocity 
attribute in the elastic half-space is calculated from vibrator ground parameters as 
discussed above. Next, a spatial smoother with 960 and 1440 meters lengths, 
respectively in the X and Y directions, was applied to the derived velocity attribute. 

15 Finally, the attribute's values were kriged in order to form an exhaustive data set. 
Step 4) VALVE model building 

A 3D near-surface velocity model was built from surface to the SRD for 
subsequent seismic statics calculation. The first model (3D kriged model) was produced 
using 3D kriging of uphole data. Uphole data were loaded to the 3D kriging engine as 
20 interval or average velocity logs. Therefore, all uphole samples contributed to the model 
regardless of their penetration depths. The second model (integrated model) was 
produced by integration of uphole velocity using collocated cokriging with the velocity 
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attribute derived from vibrator ground parameters down to 50 meters below surface and 
then combined with velocities from the first model down to the SRD (Fig. 9). The 50 
meter penetration depth of vibrator data was estimated based on numerical correlation 
between uphole and derived velocity attribute data. 

Step 5) Data and VALVE Model Consistency Checks 

Vibrator performance control data provide indications about the interaction 
between vibrator and ground. Therefore, it is prudent to investigate their repeatability 
over time. Fig. 8 shows that VP's along the north-south lines are acquired twice at 
different times. Therefore, this data set can be used to make a repeatability experiment. 
. Fig. 10 shows a cross-plot of the derived velocity estimates measured independently two 
times. • The horizontal axis of this plot represents measurements made at time 1 and the 
vertical axis represents measurements at time 2. There is a strong repeatability with a 
very high correlation coefficient. This further supports the validity of these 
measurements and their primary relations to the ground physical properties. 

It can be observed from Fig. 9 that, besides honoring the uphole data, the 
integrated velocity map also honors the features of the P-wave velocity attribute map. 
This gives the map a great deal of detail compared to a smoother map obtained using 
kriging of uphole data. 

Comparisons between die VALVE layer velocity model and the control upholes 
are shown in Fig 12. This clearly shows that the VALVE model has effectively predicted 
near surface velocities in areas with missing measurements. 

Applicatio n to Seismic Data Processing Static corrections (statics) are time shifts 
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applied to seismic traces aiming to produce a time section that is as free as possible from 
poor imaging and apparent structural features resulting from topography and near surface 
geologic variations. Statics, which primarily require near-surface velocity for calculation, 
are of vital importance for seismic reflection data processing and interpretation in land 

5 and transition zone. 

The 3D kriged and the integrated models plus a standard near-surface model were 
used in stacking the seismic data. The standard model is built based on traditional 
interpolation between uphole average velocity measurements from surface to the SRD. 
Figures 12 and 13 each show three stacks for a seismic line stacked using the previously 

10 described three models. In both figures the standard model stack shows apparent 

problems of poor continuity and structures that propagate up to the surface as compared 
to the two geostatistical models. This is attributed to the limited number of upholes in 
the vicinity of these lines. On the other hand, the integrated model shows more 
improvements than the 3D kriged model. It handles the medium and long wavelengths 

15 statics better. This is due to the fact that the integrated model has more accurate 
definition of the changing velocity boundaries because it relies on a densely sampled 
related data. 

A first conclusion was that geostatistical techniques proved to be useful in 
combining various types of data to estimate near surface velocity distributions. Honoring 
20 the spatial variability of the analyzed data allowed obtaining better velocity models as 
compared to the standard technique. The latter requires interpretation of the uphole data 
while the former does not, which in turn, reduces the turn around time required to build 
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a near-surface velocity model. 

Furthermore, using the VALVE technique, estimated P-wave velocity from 
vibrator measurements is introduced as an added attribute that can be used to improve 
determination of near surface velocity when it is integrated with uphole measurements. 
This attribute guides the area of influence of each uphole resulting in better handle of 
medium and long wavelength statics anomalies. The success obtained in the area outlined 
in this paper was followed by successful application of the VALYE technique in several 
other areas. 

The seismic image also compared favorably with that obtained using statics based 
on an independently derived near-surface velocity model obtained at a much greater 
effort. On the basis of these results, the integrated model resolved statics problems better 
than models that individually utilize uphole data through geostatistical and conventional 
computation techniques. 

A study area was selected, as shown in Fig. 16, and seismic data from a portion 
of the study area were processed using four different statics models. 

The first model constructed for this study used uphole data to construct a 15-iayer 
3D velocity model of the near surface using Ordinary Kriging from surface to Seismic 
Reference Datum (SRD). Then, travel times from surface to SRD were computed using 
this model. 

The second model incorporated the velocity attribute derived from the vibrator 
measurements in accordance with the present invention. Average velocity from surface 
to five iso-depth layers was geostatistically computed using collocated cokriging of each 
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iso-depth uphole velocity layer and the vibrator velocity attribute. Then, interval velocity 
was computed in five layers between the resulting velocity maps where each layer has a 
constant thickness of 10 meters. Following this step, travel times were computed for this 
50 meter earth's thickness using the integrated model. Travel times from 50 meters . 
5 below surface to the SRD were computed from uphole data only using 3D Ordinary 
Kriging, 

A third model was called the "Frozen Model." This model was built based on 
uphole average velocity measurements from surface to the SRD. An average velocity 
map was constructed based on uphole measurements that reach the SRD depth using 
10 conventional interpolation techniques such as least squares. Then this map was used to 
calculate statics by dividing thickness grip from surface to SRD by average velocity grid. 
This method did not utilize upholes which do not reach the SRD. It also produced 
unsatisfactory results when SRD goes above surface and when low velocity layers go 
deeper than SRD. 

15 The fourth statics model was called die "2-Layer Model. " This model relied on 

estimating the depth of the low near surface velocity based on uphole data. Any velocity 
lower than 1,800 m/s was considered to be low. Then high velocities were measured 
from uphole data to determine the amount of localized time shifts from the depth of the 
low velocity layer to the SRD. These shifts can be applied if SRD is above or below 

20 surface. However, in the study area this approach did not sufficiently resolve the statics 
problems. 

Therefore, an approach that relied on integrating the 2-Layer Model with the 
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geology of the area was utilized. This method required addition of about 400 pseudo 
upholes to the study area in order to control velocity estimation whenever needed. 

Fig. 17 shows the density of the real upholes, and Fig. 18 shows the density of 
the added pseudo upholes. As can be inferred, this method required a lot of 
interpretation and human intervention. It also required several iterations before 
converging to an acceptable model. Every iteration required model building followed by 
application of this model to stack the seismic data, which resulted in a lengthy time 
period requiring manpower and computer time. 

The SRD used in building the geostatistical models was 50 meters below Saudi 
Aramco's SRD to avoid the case when SRD goes above surface. 

Fig. 19 shows a seismic section along line CS-5. Four stacks of this line were 
produced using the previously described four models. The Frozen Model's stack seems 
to be poor compared to the other models. The 2-Layer Model produced results that are 
comparable to the results obtained from the 3D kriged model. 

The integrated model showed more improvements than the 2-Layer Model and the 
3D kriged model. It handled the medium and long wavelengths statics better. 
Automatic statics was performed using these statics models excluding the Frozen Model. 

Fig. 20 shows three seismic sections obtained from the three statics models after 
the application of automatic statics. It is apparent that better results are obtained from 
the integrated model. 

While the disclosed metiiod has been particularly shown and described witti 
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respect to the preferred embodiments, it is understood by those skilled in the art that 
various modifications in form and detail may be made therein without departing from the 
scope and spirit of the invention. Accordingly, modifications such as those suggested 
above, but not limited thereto are to be considered within the scope of the invention, 
5 which is to be determined by reference to the appended claims. 
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